if !d.name eq 'PS' then begin
    device,xsize=17,ysize=15,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end

@eu_read_bcprofl
@rhoprof
;!p.charsize=2.5
!p.multi=[0,2,1] & !p.charsize=1.0 
!x.margin=[7,0] & !y.margin=[5,3]
!P.Color =cgColor("black")
!P.Background =cgColor("white")

the1=reform(the[32,*,0])
;the0=2.0868e6 ; sdr_02
the0=3.4976e6  ; sun_hr_02
dtheta=theta[2]-theta[1]

nlev=20
thg=where(th GE -87 AND th LT 87)
restore,'velocity.sav'

rsinth=fltarr(l,m)
r2sinth=fltarr(l,m)
for k=0, l-1 do begin
  for j=0,m-1 do begin
    rsinth[k,j]=rr*r[k]*sin(theta[j])
    r2sinth[k,j]=(rr*r[k])^2*sin(theta[j])
  endfor
endfor
rhov=vv*rh2d*rsinth
rhow=ww*rh2d*r2sinth
stream1=fltarr(l,m)
stream2=fltarr(l,m)
for j=0,m-1 do begin
  for k=0,l-3 do begin
    kk=indgen(k+2)
    stream1[k,j] = int_tabulated(r[kk],rhov[kk,j],/double)
  endfor
endfor
  
for k=0,l-1 do begin
  for j=0,m-3 do begin
    jj=indgen(j+2)
    stream2[k,j] = int_tabulated(theta[jj],rhow[k,jj],/double)
  endfor
endfor

stream=smooth(stream1+stream2,2)

nls=12
lev_stream = 2.*max(abs(stream[*,thg]/rsinth[*,thg]))*indgen(nls)/float(nls-1)$
             -max(abs(stream[*,thg]/rsinth[*,thg]))

print,'levels Omega = ',minmax(omega)
omega_0=1.e9/(56.*24.*3600.)                 ;28
lev=grange(min(omega),max(omega),nlev)       ;28
title='!6a) !8T!D0!N=28d!6'
xlength=0.45
xinterval=0.05
xpos1=xinterval
xpos2=xinterval+xlength
cpos=[xpos1,0.04,xpos2,0.94]
bpos=[0.20,0.36,0.22,0.65]
contour,omega[*,thg],x[*,thg],y[*,thg],/fi,nl=64,/iso,pos=cpos,$
	title=title,lev=lev,xtitle='!8x!6',ytitle='!8y!6',xstyle=4,ystyle=4,/nodata
contour,omega[*,thg],x[*,thg],y[*,thg],/fi,nl=64,/iso,pos=cpos,/over
colorbar,range=[min(lev),max(lev)],pos=bpos,/vertical,$
         ytickformat='(F6.1)',yticks=2,ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)],$
         xaxis=0,char=1.1,title='!7X!6/2!7p!6 (nHz)',ANNOTATECOLOR=cgColor("black")
contour,omega[*,thg],x[*,thg],y[*,thg],nl=21,/iso,/over,lev=lev,color=cgColor("black")
plots,circle(0,0,0.96),thick=3
plots,circle(0,0,0.61),thick=3
plots,circle(0,0,0.718),thick=2,li=2
oplot,[0,0],[-0.96,-0.61],thick=3,li=0
oplot,[0,0],[0.61,0.96],thick=3,li=0
oplot,[0.4,0.4],[-1.,1.],thick=2,li=1
oplot,[0.6,0.6],[-1.,1.],thick=2,li=1
oplot,[0.8,0.8],[-1.,1.],thick=2,li=1
;
xpos1=xpos2+xinterval
xpos2=xpos1+xlength
; MERIDIONAL FLOW

lev_vv = 2.*max(abs(vv))*indgen(nlev)/float(nlev-1)-max(abs(vv))
lev = 2.*max(abs(vv))*indgen(nlev)/float(nlev-1)-max(abs(vv))
lev =grange(-10.,10.,nlev)
lev_vv=lev
cpos=[xpos1,0.04,xpos2,0.94]
bpos=[0.70,0.36,0.72,0.65]
jj=thg
contour,smooth(vv[*,jj],2),x[*,jj],y[*,jj],/fi,nl=nlev,/iso,pos=cpos,$
	xtitle='!8x!6',xstyle=4,ystyle=4,lev=lev
contour,smooth(stream[*,thg]/rsinth[*,thg],0),x[*,thg],y[*,thg],/over,lev=lev_stream,$
        max_value=0.,min_value=min(lev_stream),c_thick=2,c_linestyle=2,color=cgColor("black")
contour,smooth(stream[*,thg]/rsinth[*,thg],0),x[*,thg],y[*,thg],/over,lev=lev_stream,$
        max_value=max(lev_stream),min_value=0.,c_linestyle=0,c_thick=3,color=cgColor("black")
colorbar,range=[min(lev_vv),max(lev_vv)],pos=bpos,/vertical,$
         ytickformat='(F6.1)',yticks=2,ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)],$
         xaxis=0,char=1.1,title='!8u!D!7h!N!6 (m/s)',ANNOTATECOLOR=cgColor("black")
jj=where(th GE -88 AND th LT 88)
;partvelvec,smooth(ux[*,jj],2),smooth(uy[*,jj],2),x[*,jj],y[*,jj],/over,fraction=0.08,length=0.1
plots,circle(0,0,0.96),thick=3
plots,circle(0,0,0.61),thick=3
plots,circle(0,0,0.718),thick=2,li=2
oplot,[0,0],[-0.96,-0.61],thick=3,li=0
oplot,[0,0],[0.61,0.96],thick=3,li=0
print,';;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;'
;save,filename='velocity.sav',omega,vv,ww,x,y,r,th
;
!p.multi=0

end
